#include <stdio.h>
#include <math.h>
float *DirectLU(float *,int);

main() {
 int n,*n1;
 long i;
 FILE *fp;
 float *p,*p1,*x;
 n1=(int *)malloc(sizeof(int));
 if((fp=fopen("dat","rb+"))==0)
    printf("cannot open file\n");
 fseek(fp,0,0);
 fread(n1,sizeof(int),1,fp);
 n=(*n1);
 p=(float *)malloc(n*(n+1)*sizeof(float));
 p1=p;
 fread(p,sizeof(float),n*(n+1),fp);
 fclose(fp);
 x=DirectLU((float *)p1,n);
  for(i=0;i<n;i++) printf("x[%d]=%f\n",i,x[i]);
}

float *DirectLU(float *u,int n)
{
  int i,r,k,t;
  float *x,p;
  x=(float *)malloc(n*sizeof(float));
  
  for(r=0;r<n;r++) {
	for(i=r;i<n;i++)
	   *(x+i)=*(u+i*(n+1)+r);
	for(t=r;t<n;t++) {
	  for(k=0;k<r;k++)
	    *(x+t)-=*(u+t*(n+1)+k)*(*(u+k*(n+1)+r));
	}
	   
	for(i=k+1;i<n;i++)
	  if(fabs(*(x+i))>(fabs(*(x+k)))) k=i;
	if(k!=r) {
	  for(t=0;t<=n;t++) {
	    p=*(u+r*(n+1)+t);
	    *(u+r*(n+1)+t)=*(u+k*(n+1)+t);
	    *(u+k*(n+1)+t)=p;
	  }
	  p=*(x+r);
	  *(x+r)=*(x+k);
	  *(x+k)=p;
	}
	
	*(u+r*(n+1)+r)=*(x+r);
	for(i=r+1;i<=n;i++) {
	  for(k=0;k<=r-1;k++)
	    *(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i));
	}
	for(i=r+1;i<n;i++)
	  *(u+i*(n+1)+r)=*(x+i)/(*(x+r));
  }
  	  
  for(i=n-1;i>=0;i--) {
      for(r=n-1;r>=i+1;r--)
	  *(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*x[r];
	  x[i]=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i));
      }
      return x;
}

